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Abstract 

We present and discuss a variance-reduced stochastic particle method for simulating the 
relaxation-time model of the Boltzmann transport equation. The present paper focuses 
on the dilute gas case, although the method is expected to directly extend to all fields 
(carriers) for which the relaxation-time approximation is reasonable. The variance reduc- 
tion, achieved by simulating only the deviation from equilibrium, results in a significant 
computational efficiency advantage compared to traditional stochastic particle methods 
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in the limit of small deviation from equilibrium. More specifically, the proposed method 
can efficiently simulate arbitrarily small deviations from equilibrium at a computational 
cost that is independent of the deviation from equilibrium, which is in sharp contrast to 
traditional particle methods. 



The Boltzmann transport equation 
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where /(r, c, t) is the single-particle distribution function p], [^] denotes the collision op- 
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erator, a = (a x ,a y ,a z ) is the acceleration due to an external field, r = (x,y,z) is the position 
vector in physical space, c = (c x , c y , c z ) is the molecular velocity vector, and t is time, is used to 
describe (under appropriate conditions) transport processes in a wide variety of fields [2] includ- 
ing dilute gas flow p], phonon [3], electron neutron [5] and photon transport [6]. Recently, 
it has received renewed attention in connection to micro- and nano-scale science and technology 
where transport at lengthscales of the order of, or smaller than, the carrier mean free path is 
frequently considered (e.g. nanoscale solid-state heat transfer pJ|S]) 

Numerical solution of the Boltzmann equation remains a formidable task due to the com- 
plexity associated with the collision operator and the high dimensionality of the distribution 
function. Both these features have contributed to the prevalence of particle solution methods, 
which are typically able to simulate the collision operator through simple and physically intu- 
itive stochastic processes while employing importance sampling, which reduces computational 



cost and memory usage [9|. Another contributing factor to the prevalence of particle schemes 
is their natural treatment of the advection operator, which results in a numerical method that 
can easily handle and accurately capture traveling discontinuities in the distribution function 
[9]. An example of a particle method is the direct simulation Monte Carlo (DSMC) [TU] which 
has become the prevalent simulation method for dilute gas flow. 

One of the most important disadvantages of particle methods for solving the Boltzmann 
equation derives from their reliance on statistical averaging for extracting field quantities from 
particle data. When simulating processes close to equilibrium, thermal noise typically exceeds 
the available signal. When coupled with the slow convergence of statistical sampling (statistical 
error decreases with the square root of the number of samples), this often leads to computation- 
ally intractable problems [IT] . For example, to resolve a flow speed of the order of 1 m/s to 1% 
statistical uncertainty in a dilute gas, on the order of 5 x 10 8 independent samples are needed 

In a recent paper, Baker and Hadjiconstantinou have shown [pj that this rather severe limita- 
tion can be overcome with a form of variance reduction achieved by simulating only the deviation 
from equilibrium. By adopting this approach, it is possible to construct Monte Carlo simulation 
methods that can capture arbitrarily small deviations from equilibrium at a computational cost 
that is independent of the magnitude of this deviation. This is in sharp contrast to regular 
Monte Carlo methods, such as DSMC, whose computational cost for the same signal-to-noise 
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ratio increases sharply [TT] as the deviation from equilibrium decreases. 

The work in Refs [9], [12j [131 E] focused on the Boltzmann equation for hard spheres and the 
associated hard-sphere collision operator. The complexity associated with this collision operator, 
as well as others in related fields, has prompted scientists to search for simplified models; one 
particularly popular model is the relaxation-time approximation pfl [2] 

alt 

where f loc is the local equilibrium distribution function and r is a relaxation time. Despite the 
approximation involved, this collision model has enjoyed widespread application in a variety of 
disciplines concerned with transport processes [H El IH El EE EE] • 

In response to this widespread use, in the present paper, we present a variance-reduced parti- 
cle method for simulating the Boltzmann equation under the relaxation-time approximation. To 
focus the discussion, we specialize our treatment to the dilute gas case; however, we hope that 
this exposition can serve as a prototype for development of similar techniques in all fields where 
the relaxation-time approximation is applicable. Within the rarefied gas dynamics literature, 
the relaxation-time approximation is known as the BGK model pQ. 

In the interest of simplicity, in the present paper we assume r ^ r(c) and that no external 
forces are present. The first assumption can be easily relaxed, as discussed below. External 
fields also require relatively straightforward modifications to the algorithm presented below. 

As discussed in previous work [9], a variance-reduced formulation is obtained by simulating 

4 



mil 



- (/ - f° C ) 



(2) 



only the deviation f d = f — f e from an arbitrary, but judiciously chosen, underlying equilibrium 
distribution f e . In other words, computational particles represent the deviation from equilibrium 
and, as a result, they may be positive or negative, depending on the sign of the deviation from 
equilibrium at the location in phase space where they reside. As in other particle schemes [10J, 
in the interest of computational efficiency, each computational deviational particle represents an 
effective number N e g of physical deviational particles. 

A dilute gas in equilibrium is described by a Maxwell-Boltzmann distribution, leading to a 
local equilibrium distribution 

/ k ^-p("^) (3) 
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that is parametrized by the local number density n\ oc = n; oc (r, t), the local flow velocity u; oc = 
u ioc( r it), and the most probable speed q oc = ^2k^Ti oc jm based on local temperature Ti oc = 
T/ oc (r, t). Here, kj, is Boltzmann's constant and m is the molecular mass. In the work that 
follows, the underlying equilibrium distribution (f e ) will be identified with absolute equilibrium 

where uq is a reference (equilibrium) number density and Co = ^/2kbTo/m is the most probable 
molecular speed based on the reference temperature To. This choice provides a reasonable 
balance between generality, computational efficiency and simplicity. Other choices are of course 
possible and, depending on the problem, perhaps more efficient. However, care needs to be taken 



if a spatially varying or time dependent underlying equilibrium distribution is chosen since this 
results in a more complex algorithm [T3j [H]. 

Particle methods typically solve the Boltzmann equation by applying a splitting schemep in 
which molecular motion is simulated as a series of collisionless advection and collision steps of 
length At. In such a scheme, the collisionless advection step integrates 

by simply advecting particles for a timestep At, while the collision step integrates 

(6) 
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by changing the distribution by an amount [^] coll ( r ; c , t)At. Spatial discretization is introduced 
by treating collisions as spatially homogeneous within (small) computational cells of volume V ce ii- 
Our approach retains this basic structure; the particular form of these steps can be summarized 
as follows. 

Advection step: It can be easily verified that when the underlying equilibrium distribution is 
not a function of space or time, as is the case here, the advection step for deviational particles is 
identical to that of physical particles [i.e. f d is also governed by equation (jSJ) during the advection 
step]. Boundary condition implementation, however, differs somewhat because the mass flux to 
boundaries is now split into a deviational contribution and an equilibrium contribution. A more 



■""Note that a symmetrized algorithm provides higher-order accuracy-see Ref. |18) and references therein. 



extensive discussion, as well as algorithmic details can be found in 
Collision step: The variance-reduced form of ([6]) can be written as 



dt 



= ~ (f° c ~ F) - -f d (7) 
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Within each computational cell we integrate equation (EJ) using a two-part process. This integra- 
tion requires local (cell) values of various quantities, denoted here by hats, which are updated 
every timestep by sampling the instantaneous state of the gas. 

In the first part we remove a random sample of particles by deleting particles with probability 
At/f to satisfy 

f*(t + At) = f d (t) - —f d (t) (8) 

T 

In our implementation, this is achieved through an acceptance-rejection process which can also 
treat the case f = f (c). 

In the second part, we create a set of positive and negative particles (using an acceptance- 
rejection process) to satisfy 

f d (t + At) = f d (t + At) + - \f loc {t) - f] At (9) 

t L J 

This step can be achieved by the following procedure. Let c c be a (positive) value such that 
f loc (c) — F(c) is negligible for ||c||i > c c , where || • ||i is an L 1 -norm. Furthermore, let A max 
bound |/ Zoc (c) — F(c)\ from above. Then, repeat M c times: 

1. Generate uniformly distributed, random velocity vectors c with ||c||i < c c . 



2. If |j c (c) — F(c)\ > TZA max , create a particle with velocity c, at a randomly chosen 
position within the cell and sign sgn[/ ioc (c) —F(c)]. Here, 71 is a random number uniformly 
distributed on [0,1]. 

To find Af c we note that the number of particles (of all velocities and signs) that should be 
generated in a cell to obtain the proper change in the distribution function is 



cfccfr = — / |f oc -T|d 3 c (10) 
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where V ce ii is the cell volume. The (expected) total number of particles ultimately generated by 
the above algorithm is 

r Ls \f loc -F\d 3 c , . 

By equating the two expressions we obtain Ac = 8 At A max V ce ucl/(f N e o). 

We have verified the above algorithm using a variety of test cases; some representative results 
are presented below. Figure [1] shows a comparison between numerical solution of the BGK 
model of the Boltzmann equation [17J and our simulation results for the heat flux, q, between 
two parallel, infinite, fully-accommodating plates at slightly different temperatures (T and 
T\ = To + AT, AT > 0) and a distance L apart. The figure compares the heat flux normalized 
by the free-molecular (ballistic) value \qf m \ = PoC Q AT / \^/txTq) as a function of a Knudsen 
number [18] k = cqt/L, where Pq = riokbTo is the equilibrium gas pressure. The agreement 
is excellent. The simulations used approximately 50,000 particles yielding a relative statistical 
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uncertainty^ of less than 0.5%. These simulations were performed at AT = 2K, although as 
shown below, the cost is expected to be independent of the magnitude of AT in the limit of 
small deviation from equilibrium. We also performed DSMC simulations of the BGK model 
using AT = 10K and otherwise identical discretization and sampling parameters; AT = 10K 
was chosen as a compromise between best performance and a deviation from equilibrium that is 
small enough for the linearized conditions, and thus the benchmark results of [T7], to be valid. 

Figure [2] shows a comparison between a numerical solution of the linearized BGK model 
of the Boltzmann equation [19] and our simulation results for pressure-driven flow (for small 
pressure gradients). Under linear conditions, pressure-driven flow can be described [HI] by 
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where k = (1/ Pq) dP/dz is the normalized pressure gradient (here assumed to be in the 
z— direction), and x is the channel transverse direction. For k < 0, the term —kc z F can be 
included in our formulation as a source of \k\ At f >Q c z Fd 3 c = \n\ Atn c /(2^) positive and 
an equal number of negative deviational physical particles per unit volume drawn from the 
distribution 2^pn\c z \F j cq for c z > (— oo < c x , c y < oo) and c z < (— oo < c x , c y < oo), respec- 
tively. Figure [2] shows the normalized flowrate Q = —pucq/(Pqk,L) as a function of k, where p is 



2 The relative statistical uncertainty of a fluctuating hydrodynamic quantity is defined [llj as the standard 



deviation of this quantity normalized by its characteristic value 

3 To obtain the number of computational particles we divide by N e g. 
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the gas density and u is the average flow velocity (averaged across the channel width). Excellent 
agreement is observed. 

As stated above and shown in [S1H21II1], this class of deviational methods exhibit statistical 
uncertainties that scale with the local deviation from equilibrium thus allowing the simulation 
of arbitrarily low deviations from equilibrium at a cost that is independent of this deviation. 
Here we demonstrate this feature by studying the statistical uncertainty of the temperature in a 
problem involving heat transfer. Specifically, figure [3] shows the relative statistical uncertainty 
in the temperature (cry) as a function of the normalized wall temperature difference (T% — T )/T 
in the heat transfer problem discussed above, for k — 1 and T = 273K; in evaluating o"t, the 
characteristic value for temperature was taken to be the difference T\ — Tq. The standard devia- 
tion is measured from two computational cells in the middle of the computational domain, each 
containing approximately 950 particles. The figure shows that, for small T\ — T , the relative sta- 
tistical uncertainty remains independent of this quantity in sharp contrast to "non-deviational" 
methodsj. Moreover, the variance reduction achieved is such that significant computational 
savings are expected for (Ti — T )/T < 0.1. 

The algorithm described above imposes no restrictions on the magnitude of f d , although it 
is expected that the deviational approach will be significantly more efficient than traditional 



4 The statistical uncertainty of "non-deviational" methods (for the same number of particles per cell) is 
estimated using equilibrium statistical mechanics [11] . This approximation is very reasonable for small deviations 
from equilibrium |llj . 
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approaches when f d is small. If f d is sufficiently small for linearization to be appropriate, under 
some conditions, significant gains in computational efficiency can be achieved by taking the 
following into consideration. Under linear conditions, for the present model, we can write 



where u = ni oc /n — 1 and = T/ oc /T — 1. This representation can be very useful for improving 
the computational efficiency of update For example, for isothermal constant density flows, 
particles can be generated from a combination of a normal distribution and analytic inversion 
of the cumulative distribution function, which is significantly more efficient than acceptance- 
rejection. Alternatively, (fT3l) provides a means of obtaining tight bounds for \f loc — F\ and thus 
reducing the number of rejections if the acceptance-rejection route is followed. 

We have presented an efficient variance-reduced particle method for solving the Boltzmann 
equation in the relaxation-time approximation. The method combines simplicity with a num- 
ber of desirable properties associated with particle methods, such as robust capture of traveling 
discontinuities in the distribution function and efficient collision operator evaluation using impor- 
tance sampling [9], without the high relative statistical uncertainty associated with traditional 
particle methods in low-signal problems. In particular, the method presented here can capture 
arbitrarily small deviations from equilibrium for constant computational cost. More sophisti- 
cated techniques with spatially variable underlying equilibrium distribution [T3J [TJ] are expected 
to increase computational efficiency by reducing the number of deviational particles required to 
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describe the local state of the gas. One such technique is described in [2D] . 
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Figure 1: Comparison between the numerical solution of Bassanini et al. [17] and our simulation 
results (circles) for the heat flux between two parallel, infinite, fully-accommodating plates. 
Some numerical solution data for k > 1 have been transcribed from Ref. [I]. 
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Figure 2: Comparison between the numerical solution of Cercignani and Daneri [19] and our 
simulation results (circles) for the normalized flow rate in pressure-driven flow between two 
parallel, infinite, fully-accommodating plates. 
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Figure 3: The relative statistical uncertainty in temperature as a function of {T\ — T )/T for 
To = 273K. Simulation results are presented and compared to DSMC (indicated) which serves 
as a canonical case of a "non-deviational" method. The dashed line is obtained from the theory 
presented in [TT], while stars denote actual DSMC simulation results 
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